function hilbert_matrix_solver()
    % 测试不同维度
    n_values = [6, 8, 10];
    omega_values = [1, 1.25, 1.5];
    
    for n = n_values
        fprintf('\n========= n = %d 的计算结果 =========\n', n);
        
        % 生成希尔伯特矩阵和对应的b
        H = hilb(n);
        % 计算条件数
        cond_num = cond(H);
        fprintf('希尔伯特矩阵的条件数: %.2e\n', cond_num);
        
        x_true = ones(n, 1);
        b = H * x_true;
        
        % Jacobi迭代法
        fprintf('\nJacobi迭代法:\n');
        [x_jacobi, iter_jacobi, err_jacobi] = jacobi_method(H, b, n);
        
        % SOR迭代法
        fprintf('\nSOR迭代法:\n');
        for omega = omega_values
            fprintf('\nω = %.2f:\n', omega);
            [x_sor, iter_sor, err_sor] = sor_method(H, b, n, omega);
        end
    end
end

function [x, iter, err] = jacobi_method(A, b, n)
    tol = 1e-6;
    max_iter = 2000;    % 增加最大迭代次数
    x = ones(n, 1);     % 初始值改为全1向量
    omega = 0.5;        % 添加松弛因子以提高稳定性
    
    D = diag(diag(A));
    L = tril(A, -1);
    U = triu(A, 1);
    
    % 预计算迭代矩阵
    D_inv = diag(1./diag(A));
    
    for iter = 1:max_iter
        x_new = (1-omega)*x + omega*(D_inv * (b - (L + U) * x));
        err = norm(x_new - x, inf) / norm(x_new, inf);
        
        if err < tol
            x = x_new;
            break;
        end
        x = x_new;
        
        % 检查发散
        if any(isnan(x)) || any(isinf(x))
            fprintf('迭代发散\n');
            break;
        end
    end
    
    rel_err = norm(x - ones(n,1), inf) / norm(ones(n,1), inf);
    
    if iter < max_iter && ~any(isnan(x)) && ~any(isinf(x))
        fprintf('收敛于 %d 次迭代，相对误差 = %.2e\n', iter, rel_err);
    else
        fprintf('未在%d次迭代内收敛，相对误差 = %.2e\n', max_iter, rel_err);
    end
    fprintf('解向量的前5个分量: %.6f %.6f %.6f %.6f %.6f ...\n', x(1:min(5,n)));
end

function [x, iter, err] = sor_method(A, b, n, omega)
    tol = 1e-6;
    max_iter = 2000;    % 增加最大迭代次数
    x = ones(n, 1);     % 初始值改为全1向量
    
    for iter = 1:max_iter
        x_new = x;
        for i = 1:n
            sum1 = A(i,1:i-1) * x_new(1:i-1);
            sum2 = A(i,i+1:n) * x(i+1:n);
            x_new(i) = (1-omega)*x(i) + omega*(b(i) - sum1 - sum2)/A(i,i);
        end
        
        err = norm(x_new - x, inf) / norm(x_new, inf);
        
        if err < tol
            x = x_new;
            break;
        end
        x = x_new;
        
        % 检查发散
        if any(isnan(x)) || any(isinf(x))
            fprintf('迭代发散\n');
            break;
        end
    end
    
    rel_err = norm(x - ones(n,1), inf) / norm(ones(n,1), inf);
    
    if iter < max_iter && ~any(isnan(x)) && ~any(isinf(x))
        fprintf('收敛于 %d 次迭代，相对误差 = %.2e\n', iter, rel_err);
    else
        fprintf('未在%d次迭代内收敛，相对误差 = %.2e\n', max_iter, rel_err);
    end
    fprintf('解向量的前5个分量: %.6f %.6f %.6f %.6f %.6f ...\n', x(1:min(5,n)));
end